Performance Above Random Expectation

Performance Above Random Expectation (PARE)

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)
library(AUC)

## AUC 0.3.0

## Type AUCNews() to see the change log and ?AUC to get an overview.

dir.create("Figures", showWarnings = FALSE)

defaultWidth <- 8
defaultHeight <- 4.5
narrowWidth <- 4

saveFigure <- function(description1, description2, width=defaultWidth, height=defaultHeight)
{
  ggsave(paste("Figures/", gsub(" ", "_", description1), "__", description2, ".png", sep=""), width=width, height=height, units="in")
}

generateNumericScores = function(means, randomSeed=999, n=numRandomValues)
{
  if (length(means) < 1)
    stop("Must provide at least one mean.")

  set.seed(randomSeed)

  scores = NULL
  for (i in 1:length(means))
    scores = c(scores, rnorm(n / length(means), mean=means[i]))

  return(standardize(scores))
}

generateNumericScoresImbalanced = function(means, randomSeed=0, n=numRandomValues)
{
  if (length(means) != 2)
    stop("Must provide exactly two means.")
  
  set.seed(randomSeed)

  scores = rnorm(n * imbalanceProportion, mean=means[1])
  scores = c(scores, rnorm(n * (1 - imbalanceProportion), mean=means[2]))

  return(standardize(scores))
}

standardize = function(x)
{
  (x - min(x)) / (max(x) - min(x))
}

numRandomValues = 1000
imbalanceProportion = 0.95

actual = factor(c(rep(0, numRandomValues / 2), rep(1, numRandomValues / 2)))
actualImbalanced = factor(c(rep(0, numRandomValues * imbalanceProportion), rep(1, numRandomValues * (1 - imbalanceProportion))))

randomScores = generateNumericScores(0)
mediumSignalScores = generateNumericScores(c(0, 1))
signalScores = generateNumericScores(c(0, 10))
oppositeScores = generateNumericScores(c(10, 0))
aFewGoodScores = standardize(c(rnorm(990), rnorm(10, mean=10)))

mediumSignalScoresImbalanced = generateNumericScoresImbalanced(c(0, 1))
signalScoresImbalanced = generateNumericScoresImbalanced(c(0, 10))
oppositeScoresImbalanced = generateNumericScoresImbalanced(c(10, 0))

plotBarPlot <- function(values, main, width=defaultWidth)
{
  plotData <- data.frame(Scores=values)

  print(ggplot(plotData, aes(x=Scores)) + geom_bar() + theme_bw() + theme(legend.key = element_blank(), legend.title=element_blank()) + xlab("") + ylab("") + ggtitle(main))
  saveFigure(main, "Barplot", width=width)
}

plotHist <- function(values, main, width=defaultWidth)
{
  plotData <- data.frame(Scores=values)

  print(ggplot(plotData, aes(x=Scores)) + geom_histogram() + theme_bw() + theme(legend.key = element_blank(), legend.title=element_blank()) + xlab("") + ylab("") + ggtitle(main))
  saveFigure(main, "Histogram", width=width)  
}

plotBoxPlot <- function(scores, actual, main, width=defaultWidth)
{
  plotData <- data.frame(Scores=scores, Actual=actual)

  print(ggplot(plotData, aes(factor(Actual), Scores)) + geom_boxplot() + theme_bw() + theme(legend.key = element_blank(), legend.title=element_blank()) + xlab("") + ylab("") + ggtitle(main))
  saveFigure(main, "BoxPlot", width=width)
}

plotBarPlot(actual, main="Class Labels", width=narrowWidth)

plotBarPlot(actualImbalanced, main="Imbalanced Class Labels", width=narrowWidth)

plotHist(randomScores, main="Random Scores")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

plotHist(mediumSignalScores, main="Medium Signal Scores")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

plotHist(signalScores, main="Signal Scores")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

plotHist(oppositeScores, main="Opposite Scores")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

plotHist(aFewGoodScores, main="A Few Good Scores")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

plotBoxPlot(randomScores, actual, main="Random Scores")

plotBoxPlot(mediumSignalScores, actual, main="Medium Signal Scores")

plotBoxPlot(signalScores, actual, main="Signal Scores")

plotBoxPlot(oppositeScores, actual, main="Opposite Scores")

plotBoxPlot(aFewGoodScores, actual, main="A Few Good Scores", width=narrowWidth)

plotHist(mediumSignalScoresImbalanced, main="Medium Signal Scores Imbalanced")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

plotHist(signalScoresImbalanced, main="Signal Scores Imbalanced")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

plotHist(oppositeScoresImbalanced, main="Opposite Scores Imbalanced")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

plotBoxPlot(mediumSignalScoresImbalanced, actualImbalanced, main="Medium Signal Scores Imbalanced")

plotBoxPlot(signalScoresImbalanced, actualImbalanced, main="Signal Scores Imbalanced")

plotBoxPlot(oppositeScoresImbalanced, actualImbalanced, main="Opposite Scores Imbalanced")

plotROC = function(actual, scores, main, width=defaultWidth)
{
  rocResult = roc(scores, actual)
  aucValue = auc(rocResult)
  
  main2 = paste(main, "\n", "(AUC = ", format(aucValue, digits=3, nsmall=3), ")", sep="")
  
  plotData <- data.frame(TPR=rocResult$tpr, FPR=rocResult$fpr)
  print(ggplot(plotData, aes(x=FPR, y=TPR)) + geom_abline(color="darkgrey", linetype="dashed") + geom_line() + theme_bw() + theme(legend.key = element_blank(), legend.title=element_blank()) + xlab("False positive rate") + ylab("True positive rate") + ggtitle(main2))
  saveFigure(main, "ROC", width=width)
}

plotROC(actual, randomScores, "Random Scores")

plotROC(actual, mediumSignalScores, main="Medium Signal Scores")

plotROC(actual, signalScores, main="Signal Scores")

plotROC(actual, oppositeScores, main="Opposite Scores")

plotROC(actual, aFewGoodScores, main="A Few Good Scores", width=narrowWidth)

plotROC(actualImbalanced, randomScores, main="Random Scores Imbalanced")

plotROC(actualImbalanced, mediumSignalScoresImbalanced, main="Medium Signal Scores Imbalanced")

plotROC(actualImbalanced, signalScoresImbalanced, main="Signal Scores Imbalanced")

plotROC(actualImbalanced, oppositeScoresImbalanced, main="Opposite Scores Imbalanced")

calculateAccuracy <- function(labels, predictors, threshold)
{
  sum(as.integer(predictors > threshold) == labels) / length(labels)
}

calculateNumCorrect <- function(labels, predictors, threshold)
{
  sum(as.integer(predictors > threshold) == labels)
}

calculateSensitivity <- function(labels, predictors, threshold)
{
  predictorsForPositiveLabels <- predictors[which(labels==1)]
  sum(predictorsForPositiveLabels > threshold) / sum(labels == 1)
}

calculateSpecificity <- function(labels, predictors, threshold)
{
  predictorsForNegativeLabels <- predictors[which(labels==0)]
  sum(predictorsForNegativeLabels <= threshold) / sum(labels == 0)
}

calculatePPV <- function(labels, predictors, threshold)
{
  if (length(predictors > threshold) == 0)
    return(0)
  
  length(intersect(which(predictors > threshold), which(labels == 1))) / sum(predictors > threshold)
}

calculateNPV <- function(labels, predictors, threshold)
{
  if (length(predictors <= threshold) == 0)
    return(0)
  
  length(intersect(which(predictors <= threshold), which(labels == 0))) / sum(predictors <= threshold)
}

calculateFscore <- function(labels, predictors, threshold)
{
  precision <- calculatePPV(labels, predictors, threshold)
  recall <- calculateSensitivity(labels, predictors, threshold)
  
  (2 * (precision * recall)) / (precision + recall)
}

checkPareInputs <- function(labels, predictors)
{
  if (!is.factor(labels))
    stop("The labels must be a factor object.")

  if (nlevels(labels) != 2)
    stop("The labels factor object must have exactly two levels.")

  if (length(intersect(c(0, 1), levels(labels))) != 2)
    stop("The factor levels of the labels object must be c(0, 1).")
  
  if (!is.numeric(predictors))
    stop("The predictors must be numeric.")
  
  if (length(predictors) != length(labels))
    stop("The length of the predictors and labels must be the same.")
  
  if (any(predictors > 1 | predictors < 0))
    stop("Predictions must fall between 0 and 1.")
}

calculateScoresAtThresholds <- function(labels, predictors, metricCalculationFunction=calculateAccuracy, metricLabel="Accuracy", numPermutations=10)
{
  checkPareInputs(labels, predictors)
  
  # Combine the labels and predictors so we can sort them together
  combined = cbind(as.integer(labels) - 1, predictors)
  combined = combined[order(labels),]

  # Then separate them back out
  classLabels = combined[,1]
  classPredictors = combined[,2]

  thresholdStepSize = (max(classPredictors) - min(classPredictors)) / (length(classPredictors) + 1)
  thresholds = seq(min(classPredictors), max(classPredictors), thresholdStepSize)
  thresholds = thresholds[2:(length(thresholds)-1)]

  actualResultAtThresholds = calculateMetricAtThresholds(classLabels, classPredictors, thresholds, metricCalculationFunction)
  
  permutedMatrix = NULL
  for (i in 1:numPermutations)
  {
    set.seed(i)
    iPermuted = calculateMetricAtThresholds(sample(classLabels), classPredictors, thresholds, metricCalculationFunction)
    permutedMatrix <- rbind(permutedMatrix, iPermuted)
  }
  permutedResultAtThresholds <- apply(permutedMatrix, 2, mean)
  
  idealPredictors = sort(classPredictors)
  idealResultAtThresholds = calculateMetricAtThresholds(classLabels, idealPredictors, thresholds, metricCalculationFunction)

  wrongPredictors = sort(classPredictors, decreasing=TRUE)
  wrongResultAtThresholds = calculateMetricAtThresholds(classLabels, wrongPredictors, thresholds, metricCalculationFunction)
  
  whichThresholdsToKeep <- which(idealResultAtThresholds != permutedResultAtThresholds & idealResultAtThresholds != wrongResultAtThresholds & wrongResultAtThresholds != permutedResultAtThresholds)
  
  # This indicates that all predictions were the same.
  if (length(whichThresholdsToKeep) == 0) {
    pareScoreAtThresholds <- permutedResultAtThresholds
  } else {
    thresholds <- thresholds[whichThresholdsToKeep]
    actualResultAtThresholds <- actualResultAtThresholds[whichThresholdsToKeep]
    permutedResultAtThresholds <- permutedResultAtThresholds[whichThresholdsToKeep]
    idealResultAtThresholds <- idealResultAtThresholds[whichThresholdsToKeep]
    wrongResultAtThresholds <- wrongResultAtThresholds[whichThresholdsToKeep]
    
    pareScoreAtThresholds <- calculatePareScores(actualResultAtThresholds, permutedResultAtThresholds, idealResultAtThresholds, wrongResultAtThresholds)
  }
  
  thresholdData <- data.frame(ScoreType="Maximal", Score=idealResultAtThresholds, Threshold=thresholds)
  thresholdData <- rbind(thresholdData, data.frame(ScoreType="Actual", Score=actualResultAtThresholds, Threshold=thresholds))
  thresholdData <- rbind(thresholdData, data.frame(ScoreType="Permuted", Score=permutedResultAtThresholds, Threshold=thresholds))
  thresholdData <- rbind(thresholdData, data.frame(ScoreType="Minimal", Score=wrongResultAtThresholds, Threshold=thresholds))
  thresholdData <- rbind(thresholdData, data.frame(ScoreType="PARE", Score=pareScoreAtThresholds, Threshold=thresholds))
  
  return(thresholdData)
}

calculatePareScores <- function(actualResultAtThresholds, permutedResultAtThresholds, idealResultAtThresholds, wrongResultAtThresholds)
{
  scoreAtThresholds <- NULL
  for (i in 1:length(actualResultAtThresholds))
  {
    actualResult <- actualResultAtThresholds[i]
    permutedResult <- permutedResultAtThresholds[i]
    idealResult <- idealResultAtThresholds[i]
    wrongResult <- wrongResultAtThresholds[i]

    actualDiff = actualResult - permutedResult
    baselineDiff = idealResult - permutedResult

    score = 0
    if (baselineDiff != 0)
      score = actualDiff / baselineDiff

    if (score < 0)
    {
      baselineDiff = permutedResult - wrongResult
      score = actualDiff / baselineDiff
    }
    
    scoreAtThresholds <- c(scoreAtThresholds, score)
  }
  
  return(scoreAtThresholds)
}

calculatePareScore <- function(thresholdData)
{
  return(mean(filter(thresholdData, ScoreType=="PARE")$Score))
}

getMainForPlot <- function(thresholdData, main)
{
  pareScore <- calculatePareScore(thresholdData)
  
  mainScore <- paste("(Score = ", format(pareScore, digits=3, nsmall=3), ")", sep="")
  
  if (main == "") {
   main = mainScore
  } else
  {
   main = paste(main, "\n", mainScore, sep="")
  }
  
  return(main)
}

plotPareScores <- function(thresholdData, main="", width=defaultWidth, height=defaultHeight)
{
  plotData <- filter(thresholdData, ScoreType!="PARE")
  
  print(ggplot(plotData, aes(x=Threshold, y=Score, group=ScoreType, colour=ScoreType)) + geom_line() + theme_bw() + theme(legend.key = element_blank(), legend.title=element_blank()) + ylab(metricLabel) + ggtitle(getMainForPlot(thresholdData, main)))
  saveFigure(main, paste("PARE_", metricLabel, sep=""), width=width, height=height)
}

plotPareThresholdSelection <- function(thresholdData, main="", width=defaultWidth, height=defaultHeight)
{
  pareScores <- filter(thresholdData, ScoreType=="PARE")$Score
  weights <- filter(thresholdData, ScoreType=="Actual")$Score
  weightedPareScores <- pareScores * weights
  #weightedPareScores <- pareScores * weights^2

  plotData <- filter(thresholdData, ScoreType=="PARE")
  plotData$Score <- weightedPareScores

  print(ggplot(plotData, aes(x=Threshold, y=Score)) + geom_line() + theme_bw() + theme(legend.key = element_blank(), legend.title=element_blank()) + ylab("PARE") + ggtitle(getMainForPlot(thresholdData, main)) + ylim(-1, 1) + geom_hline(yintercept=0, color="red", linetype="dashed"))

  saveFigure(main, paste("ThresholdSelection_", metricLabel, sep=""), width=width, height=height)
  
  # + geom_vline(xintercept=bestThreshold, color="red", linetype="dashed")
  # + theme(axis.text.y=element_blank(), axis.ticks.y=element_blank())
}

pare <- function(labels, predictors, main="", plot=TRUE, metricCalculationFunction=calculateAccuracy, metricLabel="Accuracy", width=defaultWidth, height=defaultHeight)
{
  checkPareInputs(labels, predictors)
  thresholdData <- calculateScoresAtThresholds(labels, predictors, metricCalculationFunction, metricLabel)

  if (plot)
  {
    plotPareScores(thresholdData, main, width=width, height=height)
    plotPareThresholdSelection(thresholdData, main, width=width, height=height)
  }
  
  #return(calculatePareScore(thresholdData))
}

calculateMetricAtThresholds <- function(labels, predictors, thresholds, metricCalculationFunction)
{
  sapply(thresholds, function(x) { metricCalculationFunction(labels, predictors, x) })
}

metricCalculationFunction <- calculateAccuracy
metricLabel <- "Accuracy"
#metricCalculationFunction <- calculateSensitivity
#metricLabel <- "Sensitivity"
#metricCalculationFunction <- calculateSpecificity
#metricLabel <- "Specificity"
#metricCalculationFunction <- calculatePPV
#metricLabel <- "PPV"
#metricCalculationFunction <- calculateNPV
#metricLabel <- "NPV"
#metricCalculationFunction <- calculateFscore
#metricLabel <- "F-score"

pare(actual, randomScores, main="Random Scores", metricCalculationFunction=metricCalculationFunction, metricLabel=metricLabel)

metricCalculationFunction <- calculateAccuracy
metricLabel <- "Accuracy"
width <- defaultWidth
#metricCalculationFunction <- calculateSensitivity
#metricLabel <- "Sensitivity"
#width <- narrowWidth
#metricCalculationFunction <- calculateSpecificity
#metricLabel <- "Specificity"
#width <- narrowWidth

pare(actual, randomScores, main="Random Scores", metricCalculationFunction=metricCalculationFunction, metricLabel=metricLabel, width=narrowWidth)

pare(actual, mediumSignalScores, main="Medium Signal Scores", metricCalculationFunction=metricCalculationFunction, metricLabel=metricLabel, width=width)

pare(actual, signalScores, main="Perfect Scores", metricCalculationFunction=metricCalculationFunction, metricLabel=metricLabel, width=narrowWidth)

pare(actual, oppositeScores, main="Opposite Scores", metricCalculationFunction=metricCalculationFunction, metricLabel=metricLabel, width=width)

pare(actual, aFewGoodScores, main="A Few Good Scores", metricCalculationFunction=metricCalculationFunction, metricLabel=metricLabel, width=defaultWidth)

pare(actualImbalanced, randomScores, main="Random Scores Imbalanced", metricCalculationFunction=metricCalculationFunction, metricLabel=metricLabel, width=width)

pare(actualImbalanced, mediumSignalScoresImbalanced, main="Medium Signal Scores Imbalanced", metricCalculationFunction=metricCalculationFunction, metricLabel=metricLabel, width=narrowWidth)

pare(actualImbalanced, signalScoresImbalanced, main="Signal Scores Imbalanced", metricCalculationFunction=metricCalculationFunction, metricLabel=metricLabel, width=width)

pare(actualImbalanced, oppositeScoresImbalanced, main="Opposite Scores Imbalanced", metricCalculationFunction=metricCalculationFunction, metricLabel=metricLabel, width=width)

compareAUCvPARE = function(meanOptions, iterations=20)
{
  aucScores = NULL
  pareScores = NULL

  for (i in 1:iterations)
  {
    message(paste("Iteration ", i, "/", iterations, sep=""))
    for (j in 1:length(meanOptions))
    {
      meanOption = meanOptions[j]
      predictors = generateNumericScores(c(0, meanOption), i*j)

      aucScores = c(aucScores, auc(roc(predictors, actual)))
      pareScores = c(pareScores, pare(actual, predictors, plot=FALSE))
    }
  }

  #hist(aucScores, main=paste("AUC:", format(mean(aucScores), digits=3, nsmall=3)))
  #hist(pareScores, main=paste("PARE:", format(mean(pareScores), digits=3, nsmall=3)))
  
  plotData <- data.frame(AUC=aucScores, PARE=pareScores)
  print(ggplot(plotData, aes(x=AUC, y=PARE)) + geom_point() + theme_bw() + theme(legend.key = element_blank(), legend.title=element_blank()))
  saveFigure("AUC", "PARE")
}

#compareAUCvPARE(rnorm(100), iterations=100)